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ABSTRACT 


^Two  algorithms  are  given  for  generating  gamma  distributed  random 
variables.  The  algorithms,  which  are  valid  when  the  shape  parameter  is 
greater  than  one,  use  a uniform  majorizing  function  for  the  body  of  the 
distribution  and  exponential  majorizing  functions  for  the  tails.  The 
algorithms  are  self-contained,  requiring  only  U(0,1)  variates.  Compar- 
isons are  made  to  three  competitive  algorithms  in  terms  of  marginal 
generation  times,  initialization  time,  and  memory  requirements.  Both 
algorithms  are  faster  than  existing  methods,  for  all  values  of  the  shape 
parameter. 
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1.  INTRODUCTION 

During  the  last  few  years  many  algorithms  have  been  developed  for 
generation  of  gamma  random  variables  having  density  function 

f^(x)  » ^ exp(-x)/r(a)  O^x^”,  l<a<®. 

To  the  author's  knowledge,  these  include  Ahrens  and  Dieter  (1974), 

Atkinson  and  Pearce  (1976),  Fishman  (1973),  Fishman  (1976),  Johnk  (1974), 
Greenwood  (1974),  Marsaglia  (1977),  McGrath  and  Irving  (1973),  Tadikamalla 
(1978a,  1978b),  C.S.  Wallace  (1976),  N.D.  Wallace  (1974),  and  Whittaker 
(1974).  Most  have  been  implementations  of  the  general  acceptance/rejection 
algorithm,  with  many  using  the  modification  referred  to  as  the  "squeeze" 
technique  by  Marsaglia  (1977).  The  algorithms  developed  in  this  paper  use 
the  squeeze  technique,  which  in  the  general  case  proceeds  as  follows: 

Let  f(x)  be  the  density  function  from  which  random  variates  are 
desired  and  let  t(x)  and  b(x)  be  majorizing  and  minorizing  functions  of 
f(x),  respectively  (t(x)  ^ f(x)  for  all  x and  b(x)  £ f(x)  for  all  x) . 

Then 

CO 

1.  Generate  x having  density  r(x)  t(x)//j^  t(y)dy. 

2.  Generate  v 'u  U(0,1). 

3.  If  V £ b(x)/t(x),  deliver  x. 

4.  If  V ^ f(x)/t(x),  deliver  x.  Otherwise  go  to  step  1. 

If  t(x)  fits  f(x)  well,  if  r(x)  yields  variates  quickly,  and  if  b(x)  both 
fits  f(x)  well  and  is  quick  to  evaluate,  the  squeeze  technique  yields 
variates  quickly  even  when  f(x)  is  time  consuming  to  evaluate. 


2.  The  Algorithms 


I 
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Similar  to  the  beta  algorithms  of  Schmeiser  and  Shalaby  (1977),  the 
points  of  inflection  and  the  mode  are  central  to  this  algorithm. 

Define 

Xf  =•  X2(l  “ 1/ (x^  ~ ^2)) 

1/2 

X2  =*  Max(0,  Xj  - x^  ) 

x^  =•  ct  - 1 

. 1/2 

X,  = x_  + x_ 

4 J J 

x^  = x^(l  + l/(x^  - x^)) 

Here  x^  is  the  mode,  X2  and  x^  are  the  points  of  inflection  of  f(x) 
and  Xj^  and  x^  are  the  points  at  which  the  tangent  of  f(x)  at  X2  and  x^ 
cross  the  X axis.  If  a < 2,  there  is  no  left  point  of  inflection  and 
x^  = X2  = 0.  These  points  are  illustrated  in  Figure  A. 

Figure  A About  Here 

The  simpler,  and  slower,  of  the  two  algorithms  is  described  first. 

The  algorithm  uses  a uniform  majorizing  function  for  the  body  of  the  dis- 
tribution and  an  exponential  majorizing  function  for  the  tails.  It  is 
denoted  G2PE  since  it  is  a gamma  (G)  generator  which  requires  evaluating 
the  density  function  at  two  points  (2P)  and  uses  an  exponential  (E)  major- 
izing function  for  the  tails. 

For  simplicity  f^Cx)  is  rescaled  to 

f(x)  = expCx^  Zn(x/x^)  + x^  - x] 
to  avoid  evaluating  the  gamma  function  and  to  yield  fCx^)  “ 1. 


A 
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The  majorizing  function  is 

t(x)  - t(x^)  exp(-Xj^(x  - x^))  0 < X £ X2 

■ 1 x^  < X £ x^ 

=•  f (x^)exp(-Aj^(x  - x^))  x^  < X < « 

where  - l-Cx^/x^)  and  A^^  = l-Cx^/x^)  to  make  t(x)  tangent  to  f(x)  at  x^ 
and  x^.  Several  previous  algorithms  have  used  exponential  tails,  although 
not  as  implemented  here.  Schmeiser  (1978)  discusses  the  use  of  exponential 
majorizing  functions  for  distribution  tails  in  detail. 

The  minorizing  function  is 


b(x)  = f(x2)(x  - Xj^)/(X2  - x^) 

= fCx^)  + (1  - f(x2))(x  - - x^) 

= f(x^)  + (1  - f(x^))(x^  - x)/(x^  - x^) 
- f(x^)(x^  - x)/(x^  - x^) 

Both  functions  are  shown  in  Figure  A. 


0 < X £ X2 

X2  < X < X3 

^3  < 1 

X,  < X < 1 
4 — 


Based  on  these  functions  and  the  squeeze  technique  discussed  in 
Section  1,  algorithm  G2PE  can  be  implemented  as  follows. 

Algorithm  G2PE 

Initialization 

1/2 

1.  Set  X3  " a - 1,  D = X3  , A^  » 1,  Xj^  = X2  = f2  * O' 

If  0^X3  go  to  step  2.  Otherwise  set 

X2  * X3  - D,  A^  • 1 - X3/X2,  Xj^  “ X2  + l/^L’ 

and  f2  “ f(x2). 
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2.  Set  + D,  = 1 - x^/x^,  x^  = x^  + l/X^^- 

^4  * ^1  “ ^4  “ *2’  ^2  " ^1  ~ ^2^^’ 

P3  “ P2  ■*■  ^A^^R* 

Generation 

3.  Sample  u,  v '^^  U(0,1)  and  set  u = up^. 

If  u > p^,  go  to  step  4.  Otherwise  set  x “ X2  + u. 

If  X > x^  and  v £ f^  + (x^  - x) (1  - f^)/(x^  - x^)), 

deliver  x.  If  x < x^  and  v £ f^  + (x  - X2)(l  - ^-^1 

(x^  - x^) , deliver  x.  Otherwise  go  to  step  6. 

4.  If  u > p^,  go  to  step  5.  Otherwise  set  u = (u  - Pj^)/(P2  * . 

X >•  x^  - Jln(u)/X^.  If  X < 0,  go  to  step  3.  Otherwise  set 

V * vf^u.  If  V _<  f^Cx  - x^)/(x2  - x^),  deliver  x.  Otherwise 

go  to  step  6. 

5.  Set  u * (u  - P2)/(P2  - P2)  » “ ^4  “ ^n(u)/Xj^,  v = vf^u.  If 

V ^ ^4^’'5  ~ x)/(x^  - x^)  , deliver  x. 

6.  If  £n  V ^ x^  ^.nCx/x^)  + x^  ~ x,  deliver  x.  Otherwise  go  to  step  3. 

The  second  algorithm,  denoted  G4PE  for  reasons  analogous  to  G2PE,  is 
illustrated  in  Figure  B.  The  majorizing  function  is  uniform  over  the  body 
of  the  distribution,  triangular  over  the  shoulders,  and  exponential  in  the 
tails.  The  resulting  area  under  the  majorizing  function  is  partitioned  into 
the  ten  regions  shown.  Four  regions  have  zero  probability  of  rejection.  Of 
the  remaining  six  regions,  two  require  uniform  variates,  two  require  tri- 
angular variates  and  two  require  exponential  variates. 

The  algorithm  may  be  implemented  as  follows: 
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Algorithm  G4PE 


Initialization 

1/2 

1.  Set  “ 01  - 1,  D * , Xj^  “ X2  " ^1  " ^2  * 

If  D ^ x^,  go  to  step  2.  Otherwise  set  x^  = x^  - D, 

x^  - X2(l  - l/D),  - 1 - X3/x^,  - fCx^).  and  = f(x2). 

2.  Set  x^  = X3  + D,  X3  = x^(l  + 1/D),  Aj^  - 1 - X3/X3, 

- f(x^),  and  £3  =.  fCx^).  Set  - £2(x3  - X2) , 

P2  “ ~ *3^  P].’  ^3  “ ^1^^2  “ ^1^  ^2’ 

P4  *•  f5(x5  - x^)  + P3.  P5  = (1  - P4’ 

p^  - (1  - f^)(x^  - X3)  + p^.'p^  - (f^  - £^(*2  - X3)/2  + P3. 

Pg  “ (^4  " ^^(^5  " + p^.  pg  = -f^/X^  + pg, 

PlO  * ■*■  ^9* 

Generation 

3.  Sample  u 't  U(0,1)  and  set  u = u p^^^.  If  u > p^,  go  to  step  7. 

If  u > Pj^,  go  to  step  4.  Otherwise  deliver  x = X2  + n/f^. 

4.  If  u > p^,  go  to  step  5.  Otherwise  deliver  x = X3  + (u  - p^)/!^. 

5.  If  u > P3,  go  to  step  6.  Otherwise  deliver  x = x^  + (u  - p2)/f3. 

6.  Deliver  x = x^  + (u  - P3)/f3. 

7.  Sample  w % U(0,1).  If  u > p^,  go  to  step  8.  Otherwise 

set  X - X2  + (X3  - X2)w.  If  (u  - P^Z/CPj  - P^)  f w,  deliver  x. 

Otherwise  set  v » f^  + (u  - p^)/(x3  - x^)  and  go  to  step  13.  j 

8.  If  u > pg,  go  to  step  9.  Otherwise  set  x ■ X3  + (x^  - X3)w. 

(Pg  “ ~ Ps^  — deliver  x.  Otherwise  set 

V - £4  + (u  - P3)/(x^  - X3)  and  go  to  step  13. 

— J 

J 
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10.  Set  X = - w(x^  - x^) , V = + 2w(u  - p^)/(x^  - x^)  and  go 

to  step  13. 


11.  If  u > Pg,  go  to  step  12.  Otherwise  set  u = (pg  - u)/(Pg  - Pg) • 

X = Xj^-  (2.nu)/.\j^.  If  X ^ 0,  go  to  step  3.  If  w < (X^(x^  - x)  + l)/u, 
deliver  x.  Otherwise  set  v = w f^  u and  go  to  step  13. 

12.  Set  u = (p^Q  - u)/(p^Q  - Pg),  X = x^-  (inu)/Aj^.  If 

w < (X_(x,  - x)  + l)/u,  deliver  x.  Otherwise  set  v = w f-  u. 

K.  j J 

13.  If  £n  V > f(x),  go  to  step  3.  Otherwise  deliver  x. 
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3 . COMPUTATIONAL  RESULTS 


Based  on  the  findings  of  Cheng  (1976),  Marsaglia  (1977)  and 
Tadikamalla  (1978b),  it  appears  that  the  three  fastest  existing  algo- 
rithms are  to  be  found  in  these  three  papers.  Using  the  names  used 
in  the  above  papers,  Cheng's  algorithm  is  denoted  by  GB,  Marsaglia 's 
algorithm  is  denoted  RGAMA,  and  Tadikamalla ' s algorithm  is  denoted  by 
GAMMA. 

The  table  compares  G2PE,  G4PE,  GB,  RGAMA  and  GAMMA  in  terms  of 
generation  time  per  variate,  initialization  time,  and  memory  require- 
ments. The  times  are  based  on  the  generation  of  10,000  variates  and 
are  accurate  to  within  .02  milliseconds.  The  algorithms  were  coded  in 
FORTRAN  on  the  CDC  Cyber  72  at  Southern  Methodist  University.  The  uni- 
form variates  were  generated  by  the  relatively  fast  RANF  internal  to 
the  FTN  compiler. 


Insert  Table  About  Here 


It  is  clear  that  all  five  algorithms  are  robust  to  changes  in  a, 
with  all  algorithms  but  GAMMA  being  slightly  faster  for  larger  a.  In 
this  implementation,  the  algorithms  can  be  ranked  in  order  of  increas- 
ing marginal  times  as  GAPE,  G2PE,  RGAMA,  GB  and  GAMMA,  except  that 
GAMMA  is  faster  than  GB  for  a < 1.5.  G2PE  is  about  15%  faster,  and 
GAPE  is  30-A0%  faster,  than  the  previous  fastest  algorithm  RGAMA. 


J 

j 
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Marginal  Generation  Times  (in  Milliseconds) 


1 


and  Memory  Requirements 


Method 


a 

RGAMA^ 

GAMMA 

GB 

G2PE 

(^4PE 

1.0001 

.53 

.56 

.70 

.46 

,39 

1.2 

.53 

.61 

.67 

.45 

.33 

1.5 

.52 

.64 

.63 

.41 

.33 

2 

.52 

.70 

.62 

.42 

.33 

3 

.49 

.71 

.60 

.40 

.29 

5 

.49 

.75 

.56 

.37 

.28 

8 

.47 

.76 

.57 

.39 

.28 

20 

.46 

.78 

.54 

.40 

.28 

100 

.47 

.76 

.55 

00 

.26 

1000 

.45 

.72 

.53 

.38 

.26 

Set-up  Time 

.24 

.34 

.18 

.53-. 83^ 

1.01-1.57*’ 

Memory 

Requirements 

538 

316 

290  405 

566 

^As  implemented  using  the 
using  the  polar  method, 
and  decrease  the  memory 

normal  generator, 
add  approximately  .09 
requirements  by  224. 

For  the  implementation 
milliseconds  for  all  a 

Depends  upon  the  value  of  a.  The  lower  set-up  time  applies  when  a < 2 
and  the  higher  value  corresponds  to  a > 2.  ~ 

c 

Memory  requirements  include  necessary  routines  such  as  ALOG,  EXP  and 
SQRT. 
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Each  of  the  algorithms  requires  a one  time  initialization.  In  order 
of  increasing  set-up  time  the  algorithms  are  GB,  RGAMA,  GAMMA,  G2PE  and 
G4PE.  Since  the  algorithms  with  lower  marginal  times  tend  to  have  higher 
set-up  times,  a tradeoff  is  made  which  depends  upon  M,  the  required  num- 
ber of  variates  for  a fixed  a.  For  a £ 2,  RGAMA  is  fastest  for  M £ 3, 

G2PE  is  fastest  for  M = 4 or  5,  and  G4PE  is  fastest  for  M 6.  For  a > 2, 
RGAMA  is  fastest  for  M < 7 and  G4PE  is  fastest  for  M > 7.  For  no  combina- 
tion of  a and  M is  either  GAMMA  or  GB  the  fastest  algorithm. 

Memory  requirements  are  also  shown  in  the  table.  In  order  of  increas- 
ing memory  requirements  the  algorithms  are  GB,  GAMMA,  G2PE,  RGAMA,  and 
G4PE  which  includes  necessary  routines  such  as  ALOG,  EXP  and  SQRT.  However 
RGAMA  requires  a normal  variate  generator,  which  as  implemented  here  is 
algorithm  KR  given  by  Kinderman  and  Ramage  (1976)  with  memory  req-iirements 
of  289.  While  a normal  variate,  algorithm  requiring  less  memory  could  be 
used,  marginal  generation  times  would  Increase  for  RGAMA.  For  example, 
using  Marsaglia's  polar  method,  total  memory  requirements  for  RGAMA  were 
only  314,  but  marginal  generation  times  increased  approximately  .09  milli- 
second for  all  values  of  a.  Of  course  different  normal  generators  will 
result  in  various  tradeoffs  between  speed  and  memory. 

Ease  of  implementation  may  be  crudely  measured  in  lines  of  code  and 
additional  algorithms  needed.  In  ascending  order  of  lines  of  code  the 
algorithms  are  GB,  RGAMA,  GAMMA,  G2PE  and  G4PE.  The  only  additional  algorithm 
needed  is  the  normal  generator  used  by  RGAMA. 


; 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABU 


APPENDIX 


PROGRay  KilN  ( UrUT  .OUTp,' T . TAPF5=  InpL  T . T4PFf,  = oi  TpUT  ) 

^RU  CF  S CHl»lE  ISPR  J ^ L >_  1 .?j  I SChtu  F pn  vETnCnTST  I.N'  IVERSIT  Y 

To  CO'^PAhE  GE^ERaTIO^  BoLtINfs 

A = share  PAhAvfIE*'  if  Ti-E  GAmra  D I pT  P 1 8L  T T ON 

UIVENST  0 LLLPJ 

UATA  NAyt/?TAf')Ir*^^'*hSY,<CPNG#,^G?PFF,^G4PF*/ 

data  AS/l.u0  01fl.2M.5.?.,3.»8.*8.,70,.100.,100n,/ 

B 


N = lOOt'O 


DO  30  K=1 ’N 

- GQ  Jn  (1  

1 Call  GAyMA(ALPHA,BtTA,ISFFC,X) 
GO  To  15 

^ ^ A t i Cl_lA  Q fr 


GO  TO  lb 

3 Call  GP  (ALPHA  .RfTa*  ISEED.Y.) 

_ Go  TO  1 b_ _ 


4 Call  G2PE {alFha.BeTa, iSefD*X) 
GO  TO  15 

15  Sum  = SUM  ♦ ,« 

30  SUM2  s SI  m2  * X«x 


Time  * (5ECOND<x)-npE)*lOOO/N 
Sum  = SUM  / N 


20  WRITE  (h.ioi)  Name  (P)  ♦ALcwA,SUM«SUV|9»TTMF.*m 
101  Format  (»o^.a4»af  ic  , i6, 

1 n 


STOP 

End 
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o ">  ^ r>  o ■»  o <■»  o 


IlG»W* 


£CR(^OLT!^=  CG#^'  ( » » » '^T  A , I S ^ FP . u G A** « ) 

BPLCR  SCHHrjcpo  j(  I Y 1?»1C7P  SOUTHERN  MFThrDTST  U^IVfFSITY 
TC  GE'^tRATF  fTA^jrAPC  GAP^a  VARIATc«;  1«ING  *‘Af  AC  II  A^F  ccfr^F 
RFFtR^CS  HI'  APTKLF  IS  CD»»P.  AND  *'ATH,  WITM  AP  o L 1 C A T ICS 
VCL  3,PP.3?l-32'.  1G77 

A ■ shape  capa^FTtp 
X » GiNFPATFr  VAPTATF 

A PUST  iS  C-CFATfR  THAN  1/3  ASO  PSCC’'PFSr'  A .FT,  1 


DATA  P/1./ 

If  (p  .Fc,  A I cr  Tf*  1 

P » A 

5 » .33£3333/'CPT(A ) 

ZO  ■ l.-l.T530f l+F 

CC  » A * ZO**’  -. 5*( f-l .732351 )**2 
CL  • 3.*A-1. 

Cf  * l.-S*^ 


v=  THCf' 


REJECTI'N  PPPCFPI'PF  BEGINS  HFPE 


1 CALL  NQfML(  TFFFP»yi 
Z * S*X  + OF 

IF  (Z  .LE.  C.»  F-r  TP  1 
PGAPA  » .a*7**3 
i « -ALrC-(PANF{  IFFFP)  ) 

CO  « E ♦ .R*y**?  - PGAPA  + CC 
7 « 1.  - ZO/7 

IF  (CD  + CL*T*(1,+t*(  ,5+.3333333*TM  ,C-T.  c.)  f-C  TC  ? 
IF  (CP  + CL*ALrG(7/70)  .LT.  Q.)  GP  Tf]  1 

2 RGAPA  » RGAna*pfTA 
PcTlCN 

Ef  0 


SUBPOLTTNE  NPOf^AL  ( ^SFcr  , X ) 

GENEFATTPN  PF  n^■c  NaoNAL(0,l>  VARIATF  USING 
THE  aLGFPTTHH  GTVck  6Y  PINCERMAN  ANP  paPAGF 

IN  THE  jni-PNAL  CF  THE  AMERICAN  STATISTICAL  A^SPPIATICN  13/7* 

CCDEC  BV  PETEP  PONNEP  AND  POOIFIEP  PY  PCUCF  SCHMIFEF 
PARCH  1977  AND  Jl'NF  l<577  RESPECTIVpiV 


Data  TAIL/2.2U03EP67I66A71/ 

LU»PANF(ISEPr) 

IF  (IL'.Gt- .. ‘’P  AC7f  Ari?7987f  P ) GO  n 2 

PfTl!^^  Tpta^GL'LAP  VARIATE  E3  P£PCFNT  fF  Thf  7I*F 
Y»RANF(  ^ S?pf' ) 

X*TAJL*(  1.131  13U3'A4A1F0*LU  + Y-1.0) 


J 
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CETLBtJ 

2 I?(  LI. LY,.Q7?3K0'^<‘ ■'73606  ) RC  Tl  « 
r 

(“  tail  C3''Pt'T#^in^ 

r 

3 V»»A^’F(  ISc^n) 

W»»  A'^F  ( ! SE^r  ) 

Tl*  TAIL*T4TL /2,( 

t-ti-alcgcj) 

IF(V*V^T.GT.TI  ) C,n  rrj  3 
y«SCST{ 2.0*1 ) 

IF(Ll'.GF..'?Pf6f5A77O06<;<»C)  X»-> 
cETljPN 

4 IF(Ll.LT..Ocp72CP747q0463)  GO  TO  6 
r 

r FIA5T  II*.CAP  PFW'TTY 

C 

f V«^AMF(T5EFn) 

JSEFP'  ) 

7»V-W 

C LFT  V«  *'AX(U,V)  AMP  LET  W-f'IMCV.V) 

IF(V.6T.W)  GC  T IPO 

TEPP-V 

V«  u 

W*TEA‘P 

IOC  T«TAIL-.630B  34P0102'jq6C*V 

IF(L.LE..7''5EQ153If67tCl)  GO  TO  q 

OIFF«EX'’(-T*T*,f  ) 50662027463100-.  16002 51  ciPfP‘'f?* 

* (2.2i603“F67166471-A4S(T)) 

IFt  AB5vt2)4. 034240503750111. LE.OIFF)  GO  TO  o 
GO  TO  5 

6 IFJLI  .L7..';il?I27<»C2*’F703)  GO  TO  6 
C 

r SfCOrC  MFACLY  LTMEAE  rFNSlTY 

C 

7 V-RANF ( T«tFn  ) 

U»P4FF(  IScFO) 

Z«  V-V 

0 LFT  V«  ^AX(V,V)  AMP  LFT  G»PIN<V»w) 

IF(V.GT.h)  GO  TC  im 

TE»">»V 

V«W 

W«Ttf'P 

101  T«.4797274042??441*1 .10547366 102207C*w 
1F(  V.lE..372'*34O7fcf71790)  GO  TO  q 

0IFF.kX3(-T*T*.5) /2. 5C6 62927463 100-.ieCC251 9106 6563* 

♦ < 2.216n3'^‘=67164A71-APS(  T)  ) 

TF  (AP,5(Z)*.04O2fc44OA37312P.L£.0IFC)  GO  TO  c 
GO  TO  7 


r 

0 


THIRT  NFAPIV  ir^EAP  OENEITY 

6 V»PANF(  TSc'^P  ) 

W-PANF ( J 5cCP  I 
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z»v-v 

LtT  V«  tMP  LFT  W«PTK'(V»W) 

1F(\(.GT.*.)  GP  TP  IP? 

1F?*0«V 

V*'* 

W»TF^P 

T».  A797?7^0^.2??<.Al-.5«t  50713eC159AO*w 
IF(  V.15.  .90J==77<32AA2?ei7)  GH  TD  9 

r!PF«£  (-T*t*,r  ) /?,  KCf62P27<»fc3  ICO-.  ieCC2^1  <;K 665^^3* 

► { ?,216n3?«f71f6^7l-A?5 ( 7) ) 

IF  ( APS  ( ? »*.^‘'??77'^‘9''06eP6  .LF  .DIPP)  GO  TT  9 

GP  70  0 

X=T 

IP(Z.GE.C.O)  V»-y 

Of TIPN 

FND 


SU'5f5UUTI^P  GF  ( aUPP  A ,9FTi  » ♦ X ) 

_C HwiJCE  9 Cl-»^fc.ISPF(  ^‘bOUST  4t  1 978  q,?U.ThrPN  McT  P TcT  UNlV  E * 

c Gamma  VAhiATE  gfae^'pief.  pfeepence  chEmg < applIpi^  stati 

C ( 1 9 77  ).?(>,  1 .71 -7^  . ' 

_c Alpha  - f;HA ^ e paha ■|_Ef: 

C BETA  = scale  PAPAMt’lfcR 

r ISEED  = hAnOCM  ^LMfctF  SEen 

r X = GEnF-HATEC  GA^yp  vaciate 


Data  asavE  /-l ./ 

IE  (AlPH«  .eg.  ASAVK),  f?(«  Tn  100 

c 


C*«***lNlTIALl2AriCN 

r 

ACAV.r  - nlUjiA 


= 1.  / Sv.PT  (AtPHA*ALFHA  - 1.) 

a = Alpha  - i.3At'2S 


f *«**«GENERATiCi'l  OF  OMt  CttNMA  vARIATE  X 


100  ui  = Rai.f  ! ISEED  J 
U2  = PANF (ISEED) 

V = A » Ai  oe  till  y n .-t  1 
X = alpha  «•  ExP/V) 

Z = Ul«ljl<>u2 


IE  (R  * 2.50A077A  - A, 5*7  ,gE.  0.) 
lE  (R  ,L  r.  ALnQfZ) ) GC  tO  lOo 

2QO  X g BETAAX 

KETURiM 

End 


GO  TC  ?00 
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Sl'«PriT  Jh?  r.2P£(  «LPw/s»ECT/i,lSFco,yj 
C*****(»p(.C£  5(h»ci5FP  ULY  12,1<;7P  S^LTMeSK  •'FTHTf'T'^T  L»^UfP‘'ITY 
C TC  ttStSATC  A ^TA'JrAPD  '^AA'^'A  VAPIATP  rVPCKF^llAL  lAIl 

C SEJFCTIHn  aft  PpCTAnCLLAP  PEJ5CTICN  FOF  THE  PDCY  CF 

C C ISIS • PUTirv. 

r A » T^E  PHAFr  e\OA“FTFF  (A  ,GT.  1) 

<*  51  « Tt-t  GFNFCATEr  VALUE 


r 

DATA  ASAVc/-l./f VLL/1./ 

IF  (ALPF-A  .FO,  APAVC)  GP  TC  100 
r 

5cT-CP  REGT^*"  HF9F 
r 

ASAVE  » alpha 
>1  • 0. 
y?  • 0. 

F2«  C. 

>3  * AL  FH4  - 1 
C ■ SCRKY’) 

IF  (D  .CE.  VT)  GO  TH  10 

>2  - y3-C 

XLL  - 1.  - XB/y? 

XI  ■ XZ  ♦ 1 AXLI 

F2  « EX'*  (X3*ALrG(  »?/y?  ) + X3  - X21 
10  XA  « X3  + n 

XL9  «=  1.  - XB/XA 
X5  ■ XA  ♦ 1/XLP 

FA  . £Xr(X3*AirG(YA/V3)  + X3  - XA) 

“1  - X4-X2 

P2  . PI  - C2/XII 

F3  ■ P2  ♦ FA/vlc 
r 

r*#***VAi;  lATE  C-E''FPATTrv  opCCEDtPE  EEGIN'  HEPF 

r 

iro  t » PANMI’JFFD*?? 

V * FANMIFEFr!) 

C 

r RECIANGILAO  RFjPCTinN 

r 

IF  (U  .CT.  91)  GT  yr  200 
X » X2  + U 

IF  (X  .IT.  X?)  GT  TP  110 

IF  (V  .LT.  FA  ♦ (YA-Y)*( i-FA )/(XA-y3)  ) GC  TC  fCC 
6r  1C 

lie  IF  «V  .LT.  F2  ♦ {X-x^)9(l-F?)/(X3-X2))  Gr  in  50C 
GP  TC  ACO 
r 

LEFT  TAIL  GCHFOATTCsj 

’CC  IF  (L  .CT.  09)  Cn  TP  30C 
L « (I'-fl)  /(  P?-P1) 

X - X2  - Ainc-fU)  f XI  I 

(X  .lT.  0.)  GC  TC  IOC 

V » V ♦ F2  ♦ I 

IF  CV  .LT.  F?  * IX-vi ) /( X2-X1 ) ) GC  Tr  50O 
C-r  TC  AviO 
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r UGHT  Till  ppmfpaTT^M 

r 

3(C  L « {L-f2)  / <P3-«>?» 

> . - 4L'^r(L')  / YLP 

V . V * * I 

IF  (V  ,l7,  P**ry5-V) /(X5-V*)  ) GT  Tn  5CC 

c 

<*  FISAL  R'JtCTirK  joy 

r 

^CO  IFtAlCGfV)  .GT.  X3*4LrG(  y/>3 ) ♦ X3  - X)  GC  TG  KC 
src  X « X*R' TA 
RITLffN 
LNO 


5LPBGLTTNE  (AlOHA,  3ETA»  TSFFP*  V) 

r 

r PFLCf  SCHMFUFC  «ni’'^MFPN  mpTHCOIGT  IJMVrPSITY  f-AY  2A,  1<;7p 

C T3  GENcFaTp  fa^'ma  variates  WITH  w^ah  ALPhaapftA 

r ANC  V/BIAHfF  AfHSwBETAAEETA  USI^'G  TAHlv/  »ALt  AkS  ALFITPIThk 

r CESCRIPt"  lA'  GFNFRATIDN  CP  GAHMA  VarTaTF?  — jT 

G TC  APPFA®  TN  CAC^'. 

C VALIC  QKLY  Pro  alpha  ,GT,  1 

r 

DATA  ASAVE/-I./ 

IF  (ALPHA  .FC.  A^AVE)  GC  TH  103 

r 

r set  - 10  rPN^TAHTF 

r 

ASAVE  - ALPHA 
A ■ ALP^-A  - y, 

P * .5  ♦ .'♦SOFT  (A.AALPH4-3) 

C » A ♦ (l.*P)  / « 

r » (3-1.)  / (A*“) 

E ■ EXP(-A/P)  / ?. 

r 

r GEHER/TIGH  rip  rvr  VAFIAT£ 

r 

U-vi  L - F ♦ RAWFfTSfPf')  * (l.-E) 

IF  (L  .CT.  . p)  GO  TC  200 
X ■ A ♦ PWALrG {(♦('» 

IF  (X  .LT.  0.)  GO  TC  ICO 
T » A - X 
GC  TC  3tC 

?ro  X « A - pwAL nG( ?, -I  -i,n 

Y ■ V - A 

irc  L • RAMfdPPFD 

IF  (ALOC(U)  .C-T.  (A*ALCG(0*X)-X  + (Y/3)Hf ))  Cr  TC  ICC 

X ■ X*3FTA 

RFTLIRN 

EMO 


™is  PA®  K bet  WALin 


C 

SUBPOUTIaE  G^aOE  (ALi-t-A,RcTA,ic;PEC.Xi 

GHUCE  fO^'EISPW  Jt.L>  10-74  ?OLT^•FCA,  mctpO''IST  MkIvEPSIT 

c 

c 

c 

GENEPAT/CN  CF  0^  h p ''C  UCn_»AH(JOP  VARtATE  L'SimG 

TmE  FOI’m  point  MtTPCL  IaItN  FXPCnEmtiAl  T^IlS 

from  !(•,£  C'mPNA  ntNilTY  FiMCTIon 

c 

c 

r 

Alpha  = shape  pama^’eter 
dETA  = scale  oAndA/tlER 

IBEEO  = PAI.CCb;^  Mjwtr.p  cFpn 

c 

c 

_c_ 

X s GEN'EhATEC  Gan^A  \,APIa''’E 

REFEPE^CE  H.  SchmEI-SEP  SGFEZE  r'FTMCOS  FQc  C-Ai-'^A  VAPIAT' 
GENFRATIUNf  ORE'X  TtOP  rFPdPT  7annq  Jill  Y IPYR 

c_ 

Data  ajav/E  /-I  . / I XLL/-1 , / 

IF  (ALPHA  .EC.  aEaVE)  Gr  ■’’o  100 

CflflflflfliiMjT^AL]  jaTICN 

ASAVE  = ALPHA 

1. 

X1=x?sF1=fc=0, 

X3  = Alpha  - i . 

0 = -SOPT  (X 3)  . , 

IF  (D  .GE.  x3)  an  re  lo 
^2  = X3  - 0 

XI  = x?*- ( 1 .-1  */n) 

XlL  = ] .-XJ/Xl 

Fl  = EXP  (X3«ALnG  (Xi/x3)  ♦ x3  - Xi) 

F?  = rxt>  (A.nifti  PL- (X>Ayx-i^ *jrL  - x?,  _ 

2. 

10  X4  = X3  ♦ u 

xs  = xat>  { 1 . ♦!  ,/n ) 

XLR  s 1 . - X3/XR 

C 

F4  = FXP  (X3flALn6 (Xh/ X3 ) ♦ X3  - X4) 
r5  = EXP  (X3«ALne (Xb/X3)  * x3  - XR) 

c 

r 

Calculate  prcrapility  foc  Each  of  the  ten  regicns 

Pi  = F?fl  (X3-X?) L 

P2  = r4o(X'+-x'’)  ♦ Pi 

P3  = E]  T*  (X2-X1  ) ♦ Pc 

P4  = Ec«(Xb-X4l  •»  p-^ 

PS  = ( 1 ,-F2) fl (X3-x2)  ♦Pa  * 

P6  = ( ] ,-F4)  « (Xa-x3  ) 4.  pq 

P7  = tP?-Fl ) fl f X5-xi ) fl  .c  , Pa 

Pa  = (F4-F5) « fXq-Xfl ) fl ,5  ♦ P7 

P9  s -PWXUL  ♦ Ptt 

^ — E44I.-j;.FS/y.LR..A-ai» 

(;TH*<»«<AGENEPATt  ONE  Gam^a  VAPIatE  X 


quality  mcricABLi 

COl  Y fUiittl5ii£u  TO  JiLC 


GO  TO  ]'»UC 

IF  lu _*£)li, J=JLl G- 0 U'  

X = XI  * <U-F?)  / Pi 
Go  To  I'tOO 

X = Xd  ■>  (U-F1)  / F’? 

Go  TO  HJO 

IHF  TWr  N FulON'S  LS  ,Ui;_c  FfTuhi  0^  LA^H  F^FrTT< 
w = RAKKdSEEn) 

IF  (U  . GT.  P5l e 0 _ TjCL_P.tlQ 

X s X2  (a3-x2^  o a 

IF  ( (U-F'> ) / (F'=i-P‘* ) .LF.  W)  Go  TO  lAOO 

V = Fi;  ♦ IG-FA-)  / (A3-X?) 

GO  TC  130D 

IF  (b  .'’T.  P6l  Go  To  700 

X = X3  ♦ (X4-x3>  o v» 

IF  ( (P6-U ) / (F6-p=  ) .QE.  >V)  Go  TO  I4OO 

V = Fd  * (U-F5) / < Xd-«3) 

on  TO  1300 

The  TrtO  TRTAKGULflF  CErslCNS 


IF  (u  .GT.  Pfi)  GO  r,c  <;oo 
V.2  = RAWFdSEPG^ 

? 


U .or.  P7 ) Go  “0  Soo 
X)  ♦ (X2-X1>«W 


IF  'V  .Lt,  F2«W)  eo  TO  idOO 
GO  TO  )300 


V = F5  * 2*»X<»  {u"P7  ) / ( X5,  <4  ) 
GO  TO  1300 


THE  TwO  EXPCNiENUdL  RpOlONS 


(pq-U)/ (pq-pa) 

XI  - ALOG (U) /XLU 
f*  _iF-  (in 


IF  (W  .Lf.  (XLL« tXl-X ) +1 . ) /U)  GO  To  1400 

V = W»Fl«U 


1000  U = (Pi  0-U) / (P] o-P^ ' 
X = X5  - AlOG(U)/xLh 


WopS<*U 


1300  IF  (ALOGlV)  .gT,  x3* a lC G ( X/ X3 ) ♦ X3 


50  To  100 
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